Astronomy & Astrophysics manuscript no. 7132 


© ESO 2008 


February 4, 2008 





A new equation for the mid-plane potential of power law disks 
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ABSTRACT 



Aims. We show that the gravitational potential if/ in the plane of an axisymmetrical flat disk where the surface density varies as a 
power s of the radius R obeys an inhomogeneous first-order Ordinary Differential Equation (ODE) solvable by standard techniques. 
Methods. The exact derivative of the midplane potentiel in its integral form is found to be algebrically linked to the potential itself. 
Results. The ODE reads 

— - (1 + s)- = A(R), 
dR R 

where A is fully analytical. The potential being exactly known at the origin R = for any index s (and at infinity as well), the 
search for solutions consists of a Two-point Boundary Value Problem (TBVP) with Dirichlet conditions. The computating time is 
then linear with the number of grid points, instead of quadratic from direct summation methods. Complex mass distributions which 
can be decomposed into a mixture of power law surface density profiles are easily accessible through the superposition principle. 
Conclusions. This ODE definitively takes the place of the untractable bidimensional Poisson equation for planar calculations. It opens 
new horizons to investigate various aspects related to self-gravity in astrophysical disks (force calculations, stability analysis, etc.). 
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1. Introduction 

The Poisson equation plays a major role in a large variety of 
astrophysical problems where the dynamics of particles and flu- 
ids is influenced by gravitation. Like most Partial Differential 
Equations (PDEs), it is difficult to solve accurately in space and 
particularly inside matter due to its tridimensional nature and be- 
cause it needs proper boundary conditions which are not easily 
accessible in the vicinity of systems. There are many situations 
where potential values or forces are required in a single plane 
where matter is gathered through a flat disk. In such cases, the 
Poisson equation in cylindrical coordinates (R, <f>, z) 

(R-'dnRdR + R~ 2 d 2 w + 4) if/ = AnGp(R, cf>, z), 

where p is the mass density and G is the gravitation constant, 
is of no use without i) the extension of the computational grid 
along the z-direction perpendicular to the disk plane and ii) the 
additional knowledge of the gravity vertical gradient in the plane 
of the disk — the most critical point. Unfortunately, the quantity 
d 2 ; iff is not available in a simple manner. This is the reason why 
ifr is always determined through direct summation techniques. 

In this letter, we report the discovery of a new local equation 
satisfyied by the midplane potential ift in flat axially symmetri- 

XOO 
^pdz varies as a 

power law of the cylindrical radius R. As shown, it is an inhom- 
geneous first-order Ordinary Differential Equation (ODE) sub- 
ject to perfectly known Dirichlet boundary conditions. In some 
sense, this new equation circumvents the d^-problem quoted 
above by converting the genuine PDE into an ODE for planar 
calculations. Also, it is solvable by standard algorithms with 



very low computational cost scaling linearly with the number 
N of grid points (instead of A^ 2 with summation meth ods). As 
powe r laws are often met in astro physical disks (e.g. Pringle, 
Il98lb lAndrews & Williams!. I2005I) . this ODE seems very well 
suited to investigate all these cases where disk gravity is to 
be acounted for. More generally, this represents a significative 
progress in potential theory with applications in other fields of 
physics. 

2. Basic theoretical considerations 

The potential if/ due to a flat axially symmetric d isk in its plane 
at a cylindrical distance R from the center is (e.g. lDuran dl fl95l 



ifr(R) 



-2G 



f 



' —~L(a)mK(m)da, (1) 
R 



where K is the complete elliptic integral of the first kind, m is 
the modulus with 



2 ^ 

m = and < m < 1, 

a +R 



(2) 



flin > is the inner edge and a out > a m is the outer edge. We 
consider surface density profiles in the disk of the form 



1(a) = 2 0Ut | 

i a out 



(3) 
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where S out is the outer edge value and s is a real exponent. This 
choice is natural. Power laws with negative expo nents are com- 
monly met in disk theories (e.g. iPringlei [ 1981) and observa- 
tions as well (e.g. lAndrews & WilliamsL 120051) . These can also 
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be combined to describe other kinds of mass distributions (see where 
Sect. [8}. The total disk mass is then 



M d = In 



S left(c7) = - 



I,(a)ada = 27rE out a out £(A), 



wher43 



1 - A 



2+.S 



2 + 5 



(4) 



(5) 



7ix(A)m 



(14) 



A key quantity is the potential at the center i//(Q) = \f/ c which 
can easily be deduced from Eq.(fl~|i since K(0) = |. We have 



We see that the normalized potential obeys an inhomogeneous 
first-order Ordinary Differential Equation (ODE), with S\ e ft as 
second member and s and A as parameters. 

For/? > a ut (i e. outside to the disk), we proceed in a similar 
manner with u — a/R < 1. Using Eq.©, Eq.(fT]) becomes 



ifr(R) = -4GR 



if/ c — -2nG 
where (see noteQ} 



I S(o)K(m)m£/m, 

Ju ia 



(15) 



~L{a)da = -2^GS out a out ^(A), 



x(A) 



1 - A 1+ 
1 + s 



(6) 



(7) 



where u m — a m /R and u out - a oal /R. Taking the exact derivative 
of \$i and using normalized quantities, we find 



# (7/ 

-j- =(l+5)-+5 right (C7), 

dm m 



and A = ai n /a ou t is the axis ratio. Note that cases s — {-2, -1} are where 
not compatible with a disk filling the center (i.e. with a m = 0). 



5 right (ct) = _ 



nx(A)m 2 



K - -K - A 



.v+2 



(16) 



(17) 



3. Towards the ODE 

This is the same ODE as Eq.([T3l>. but with a different second 
In order to obtain the ODE, we divide the /?-axis into 3 domains, member 



For R < a m (i.e. inside the central hole), we set v = R/a < 1 and 



Finally, for a m < R < a out (i.e. within the disk), we use both 



change the modulus of the complete ejliptic integral according the variable v an d u. The potential reads 
to the identity (IGradshteyn & RyzhikL 119651) 



K 



2 V* 
1 + x 



(1 +x)K(x), 



where x < 1 . 



(8) 



-4GflE out 

out 



Since E(a) = E out v out v s , Eq.(Q} also reads 



f K(u)u l+S du 



(18) 



dv 



ifr(K) = 4G/?X 0Ut v out f °" K(y) , ( 9 ) 



and so 



— = (1 + s)- + 4G£ 
dR R 



K(Vout) 



- K(M in )M. 



(19) 



where Vi n = R/a m and v out = R/a out . The potential depends on 
five quantities: four parameters a m , a out , s and X out and a single ._. 
space variable R. The radial gradient of the potential (i.e. the Multiplying this expression by a out /^ c and using Eq.© again to 
opposite of the acceleration) is then given by an exact derivative eliminate ^out, we obtain the expression 



dtlr 4* .. d 

— — — (I + s)— +4GRH oai v* ,— 
dR 'R out out dR 



{£ 



K(v) 



dv 



(10) 



dm m 



,(m). 



(20) 



It is not necessary to perform the integration, as using 



We recognize the same ODE as above, but with yet another sec- 
ond member, namely 



d r yi{x) 

dx J yiW 



/■<vx/v= /■<3*)^-/(yi)^r. (in 

ax ax 



* inside 



(m) 



7rx(A)m 2 _ 



mK(m)-K\ — \A 



s+2 



(21) 



we directly get 



^=(1 + ^ 
dR R 

+ 4G/ffi out vf 



(12) 



K(vout) dv out K(v in ) dv u 



2 +s 
out 



,,2+s 



dR 



where dv/dR = 1/a. Let tfr = ^/^ c be the normalized potential 
and or = R/a out the normalized radial coordinate (it is equivalent 
to v out ). Using Eq.© to eliminate E out , Eq.dTZb takes the form 



4. Asymptotic properties 

It is interesting to verify that the differential equation possesses 
the right properties both at small and at large distances. A second 
order expa nsion of the term Si e f t (m) as ml A — > for a m > 
shows that dGradshtevn & Rvzhiktll965l) 

(l-A- 1 ) 

lmi^ M ( W )=-(l + ,)- — 



m\ 



(22) 



ddr dr 
dm m 



(13) 



Note that lim„ 



- In A for any A ± 0. 



As a consequence, we get the approximation 

d®-l) - m 2 (l-A°- 1 ) 



(23) 
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which admits a trivial solution, namely 
for) 



m 2 /l-A*- 1 
1 + — T 



%(A) 



1 



(24) 



As expected, the potential is quadratic with m close to the cen- 
ter, implying the linearity of the gravitational acceleration with 
m. This approximate expression is fully compatible with Eq.© 
provided the complete elliptic integral is preliminarily second- 
order expanded over the modulus as R — > 0. Note that the ODE 
is not defined at the center. This is however not a big problem 
since iff — 1 and dif/ldm = at nr — by symmetry. 

Conversely, for R » a out , the potential is expected to merge 
with that of a central point mass with mass M A . We have m — > 
as tu — > oo, hence ^(A)cr 2 5 1 -ight(TO') — > A I+2 - 1. From Eqs.(0]l 
and ([§}, we then find 



difi 
dR 



( GM A \ 



GMd 



(25) 



which must be fulfilled for any s. This has a double implication: 
ip -> -^fi and -diff/dR -> as /? oo. This is correct. 

5. The ODE in compact form 

As shown, at any point in the disk plane, iff satisfies the ODE 

^--(\ + s)^=S{^), (26) 
dxn zu 

where the second member is a piece-wise defined function 



S(m) 



S ieft(w), for < zu < A, 

Sinside(w), for A < ZU < 1, 

Sright('W), for BT > 1. 



(27) 



These three expressions are apparently different. However, 
they can be converted back into a single function by changing 
the modulus of the function K using Eq.©. From Eq.(f2]l, we 
have 



A + zu 

and so we get the expression 



and 



S(zu) 



K(m in )m in A' s 



2 yfnr 
1 + zu' 



K(m out )m out 



(28) 



(29) 



valid for any zu e [0, oo[. One thus recognizes the expres- 
sion given in the abstract, with A(zcr) = \ft c S (zu) / Wi. This ODE 
can obviously be rewritten in different equivalent forms. With 
caution, it can even be deduced by direct calculus (see the 
Appendix ). Note that E q.d26b is compatible with the famous so- 
lution of MesteJ d 1963b where the square of the disk rotational 
velocity -d\ nTn \p is uniform for s — -1 provided the disk has 



infinite mass/extension (i.e. a m = and a 



3). 



6. Gravitational forces and zero-gravity point 

Eq.j26T> can also be regarded as an algebraic relationship be- 
tween the potential iff and the radial acceleration of gravity 
gjt - -difr/dR. This is particularly attractive from a dynamical 
point of view since the force Fr = figR undergone by a particle 



with mass yu moving inside the disk plane is directly accessible 
from the potential (instead of its gradient), by the formula 



F R (R) = -\i 



t/r c 



(l + s), 

ipivT) — S (tu) 



(30) 



There is always one point inside any disk where the potential 
goes through a maximum (and the total force vanishes). Let 
■CTeq 6 [A, 1] be this equilibrium radius in units of a out ; we have 



(1 + i)(Aeq = — nT eq 5(ZiTeq)- 



(31) 



7. The ODE, in practice 



As it is well known, ODEs are easier to solve than PDEs. Here, 
tfr is exactly known in two places in the disk plane: at the cen- 
ter where iff - 1 by construction, and at infinity where iff — 
by nature. Thus, the problem is intrinsically a two-point bound- 
ary value problem (TBVP) with exactly known Dirichlet condi- 
tions. The presence of the second member S presents no special 
difficulty since special functions can be determined at computer 
precision from numerical libraries. Obviously, any finite size do- 
main [cTieft > 0, CT r ight < 00 ] can a ls° be considered, provided the 
normalized potential is accurately determined at the left bound 
lAieft = <A( CT ieft) an d right bound fright = 'A(c'right) of the interval. 
For instance, on a regular ?zr-grid made of + 1 mesh points, 
each with coordinate -nr, = m\ e f t + it where I is the mesh size 
(i.e. Ni = Wright - wi e ft) and i = 0, . . . , N, then the ODE can be 
discretized as follows 



fon - 2(2 + s)£— - h-i 



for i = l,N-l, (32) 



at second order, with iffQ = i^i e f t and tfr^ = bright - These N— 1 alge- 
braic expressions can be put in matrix form, and rapidly solved 
for iffj by standard techniques (e.g. with the Thomas algorithm). 

If a potential value ifri^ is known at a "starting" radius zD-, e ft 
(possibly the center), solving the ODE amounts to an Initial 
Value Problem (I VP) since Eq.([26*]i also writes 



iJf(V7) = l^left + 



f 



(i + a) fc=2 +W ) 



dm' . 



(33) 



A lst-order implicit scheme (slightly better than an explicit 
scheme) appropriate for s < would then give at the value 



07+1 = - 



)S M ] 
SVJi+i - (1 + s)vji 



(34) 



and so the solution can be carefully propagated up to Wright. 

We immediately see two major advantages of solving the 
ODE with these methods: first, the computing time is expected 
to be proportional to the number of grid points (instead of 
Af 2 with direct summation methods) and second, the method si- 
multaneously supplies the potential and the acceleration. Let us 
give a simple example. Figure QJ displays the potential iffi(vjj) 
when the solution of the IVP is propagated from the center 
(v7o,if/o) = (0,1) through a disk using Eq.(f34b with the fol- 
lowing parameters: Af = 1000 and vtn = 1.1 for the grid, 
and fl; n = 0.1, fl ou t = 1 (axis ratio A = 0.1) and s = -1.5 
for the disk. Figure QJ) displays the corresponding acceleration 
g m = -difffdnr obtained from Eq.(F26]i. We have compared these 
data with the 1 6-digit reference value s obtained from the split- 
ting method bv lHure & Pierensld2005l) . The relative errors are of 
the order 3 x 10~ 4 on average, which is already remarkably low. 



the mid-plane potential of power law disks 
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reference value 
— solution of the ODE 
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0,5 1 
normalized radius BJ=R/a„,„ 



__. i I i i i i L 
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Fig. 1. Normalized potential if/ obtained by solving the ODE 
(left), and acceleration of gravity deduced from Eq.(l26l> (right) 
(see text for the disk parameters and set up). Reference va l ues ar e 
computed from the splitting method bv lHure & Pierensl 12005). 
The relative errors, less than 0.1%, are not visible on the graphs. 



properties of the ODE have been discussed and a simple example 
of numerical integration through one of the most basic schemes 
has been given. As quoted, the combination of power law pro- 
files make the description of a wide variety of matter distribu- 
tions possible. This work can be analyzed in much more detail 
(like the numerical implementation of the ODE) and extended in 
several ways. For instance, it would be very interesting to con- 
sider other density profiles, the off -plane case as well as the full 
planar case by releasing the 0-invariance. These questions will 
be tackled in forthcoming papers. 
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Better results can be obtained with other schemes. A crude com- 
parison of computing times shows that, at a given precision level, 
potential values and accelerations are determined together much 
more rapidly from the ODE (by a factor ~ 40 in the present 
example) than from direct summation. This is far from being 
negligeable, especially if the potential is to be determined many 
times. 



8. Complex distributions 

The linearity of the Poisson equation enables the description of 
complex distributions of matter from individual systems through 
the superposition principle. This is particularly attractive since 
power laws represent probably the most natural (if not the sim- 
plest) set of basis functions (in particular, a basis for polynomi- 
als). So, a composite system made up of / co-planar and co-axial 
disks, each with inner edge a^j, outer edge a out j, axis ratio Ay, 
outer surface density "L ont j and power index Sj, generates at dis- 
tance R from the center a total potential 



lAtotW = ^^/(cT^iAc,, 



(35) 



where zjj = R/a 0UtJ and i/r C; = -27rG£ out;i a out>jr ^(A ; ) according 
to Eq.©. The problem then involves J independant ODEs which 
can be rapidly solved, possibly in parallel. This is illustrated in 
Fig. [TJ which shows the potential (same disk as before) for an 
exponential profile £(a) oc exp(-a/a out ) reconstructed from its 
Taylor expansion (a series of power laws). Here, i[r tol has been 
determined from the TBVP and space mapping, with J = 20 
terms (sufficient to perfectly match the exponential in the range 
[0, 1] with 16 digits). 

9. Concluding remarks 

We have demonstrated that the classical Poisson equation can 
advantageously be replaced by an ODE when potential values 
are required in the plane of an axially symmetrical flat finite size 
disk where the surface density is a power law of the radius. Some 



Appendix A: Direct derivation of the ODE 

Equation ( f26b can be derived directly from Eq.(Q]) by integrating 
the kernel over the modulus m. By reversing Eq.(f2|, we find 



a 
R 



1 + m' 



= u(m), 



(A.l) 



where m' = - m 2 is the complementary modulus. The sur- 
face density 2(a) is then an explicit function of m and R with 



2(a) = Y.(m,R) = Y, out u s (m)u ^, 
as well as the exact derivative 



da 
dm 



2Ru(m) 



1 + u(m) 



1 - u(m) 



(A.2) 



(A.3) 



Under these conditions, the midplane potential writes 
ifr(R) = -2GX out RuZ 



M s+ 2(m)|— — \K(m)mdm. (A.4) 
R dm J 



Note that the term 



R dm 



(and subsequently the whole inte- 



grand) is a function of m only. As Ru Q * t oc R l+ \ we have 



dR 'R 



(A.5) 



2G2 out 7?M out - 



± 1 r 

Using identity fJSJl, we find 



5 (m) I I K.(m)mdm 

J \Rdm' 



— = (1 + -2G2 0Ut « 0Ut « out - K( moot )m oat — 



s+h (da\ dm m 

UL K(mm)min ^" 



(A.6) 
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Since 



/ dm\ 


m 


u(m) — 1 


\dR) a ~~ 


' 2R 


u(m) + 1 



the above expression can be rearranged into 

dR y 'R 

+ 2GS out M^ ut [K(m out )m out - A' 5+I K(m in )mi n J . 

which is equivalent to Eqs.d26*i> and i 



(A.7) 

(A.8) 
(A.9) 



